# Dependency: sim.results in "sim_results_miss2.Rdata"

library(reshape)
library(ggplot2)

# Load simulation output

load("sim_results_miss2.Rdata")

J <- 1000
K<- 19

# Bias formula

bias <- function(v,true) {
  sum(v - true)/length(v)
}

true.effect <- NULL
sim.results.long <- NULL
sim.results.mtx <- NULL
est.bias <- NULL

for (k in 1:K) {

true.effect[[k]] <- sim.results.miss[[k]][1]
sim.results.long[[k]] <- sim.results.miss[[k]][-1]

sim.results.mtx[[k]] <- matrix(NA, nrow = length(sim.results.long[[k]]), ncol = J)

for (i in 1:length(sim.results.long[[k]])) { # ITERATE OVER SIMULATIONS
  
  sim.results.mtx[[k]][i,] <- as.matrix(sim.results.long[[k]][[i]][,1])
  
}

est.bias[[k]] <- apply(sim.results.mtx[[k]], 1, bias, true = true.effect[[k]][[1]][[2]])

names(est.bias[[k]]) <- names(sim.results.long[[k]])

}

# Plot bias by value of heterogeneous effect of G

est.bias.matrix <- matrix(unlist(est.bias), byrow = T, ncol = length(sim.results.long[[1]]))

colnames(est.bias.matrix) <- names(sim.results.long[[1]])

PropG_sample.vector <- c(0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95)

dev.vector <- PropG_sample.vector/0.5

rownames(est.bias.matrix) <- dev.vector

weighted.methods <- c("vATT.wweight", "vATT.wweight.reg", "vATT.wnn", "vATT.wnn.reg", "vATT.ws", "vATT.ws.reg")

data.plot <- melt(est.bias.matrix[,colnames(est.bias.matrix) %in% weighted.methods])

names(data.plot) <- c("Beta_G", "method", "bias")

# Figure 2, Panel B

ggplot(data=data.plot, aes(x=Beta_G, y=abs(bias), group=method)) +
  geom_line() + labs( x = expression(Ratio(Prop(S[2],sample)/Prop(S[2],population)))) + expand_limits(y=c(0,0.4)) + theme_bw(base_size = 18)


